function [Rho,Vp,Vs]=Dvorkin(Phi,PhiC,Clay,Sh,Pe,n)
%% Input parameters
% Phi = [0.6:-0.001:0.1]';
% PhiC = 0.4;
% Clay = 0.6;
% Sh = rand(size(Phi,1),1);
% Pe = 0.0003; % GPa

%% Define elastic moduli (GPa) of quartz, clay, hydrate, water (after Dvorkin,1999)
Kqua = 36.6; 
Gqua = 45;
Kclay = 21;
Gclay = 7;
Khyd = 7.7;
Ghyd = 3.2;
Kw = 2.25;
Gw = 0;
% Densities
Rw = 1.035;
Rhoh = 0.91;
Rhoqua = 2.66; 
Rhoclay = 2.58; 

%% Mass balances
% New porosity
Phi0 = Phi; % Original porosity
Phi = Phi0 .* (1 - Sh); % porosity after hydrate fills part of it
vqua = ((1 - Clay) .* (1 - Phi0)); % Mass of quartz (per unit of rock)
vclay = Clay .* (1 - Phi0); % Mass of clay (per unit of rock)
vhyd = Phi0 .* Sh;  % Volumetric concentration of hydrate in the rock (Helgreud et al. 1999)
% Define fraction of quartz, clay, hydrate
fqua = vqua./(1-Phi);
fclay = vclay./(1-Phi);
fhyd = vhyd./(1-Phi); 

%% Densities and pressure
% Solid-phase density
Rhom = fqua.*Rhoqua+fclay.*Rhoclay+fhyd.*Rhoh;
% Bulk density 
RhoB = (1-Phi).*Rhom + Phi.*Rw;
% % D = ([1:2:602]')./100; % depth bsf in centimeters
% % Pe = (RhoB - Rw) .* 981 .* D; % effective pressure in Pa
% % Pe = Pe .* 10.^(-9); % convert to GPa

%% Moduli
% Compute elastic moduli of solid and fluid (Hill's average solid frame)
Kf = Kw; % load bearing model = fluid is only water
Gf = 0;
Ks = 0.5.*(fqua.*Kqua + fclay.*Kclay + fhyd.*Khyd + ...
    1./(fqua./Kqua + fclay./Kclay + fhyd./Khyd));
Gs = 0.5.*(fqua.*Gqua + fclay.*Gclay + fhyd.*Ghyd + ...
    1./(fqua./Gqua + fclay./Gclay + fhyd./Ghyd));

% n = 8; % Coordination number
nu = 0.5 .* (Ks - (2/3).*Gs) ./ (Ks + (1/3).*Gs); % ??

% Hertz-Mindlin (1949)
Khm = ((Pe .* n.*n .* Gs.*Gs .* (1-PhiC).^2)./ (18.*pi.*pi .* (1-nu).^2)).^(1/3);
f1 = (5 - 4.*nu)./(10 - 5.*nu);
n2 = 3.*n.*n.*Gs.*Gs.*(1-PhiC).*(1-PhiC);
d2 = 2.*pi.*pi.*(1-nu).*(1-nu);
Ghm = f1 .* (n2 .* Pe ./ d2).^(1/3);
%% Moduli dry rock
for i=1:size(Phi,1)
    if Phi(i)<PhiC
        % Phi < PhiC
        Kdry = 1./((Phi./PhiC) ./ (Khm + (4/3).*Ghm) + ((1 - Phi./PhiC)./(Ks + (4/3).*Ghm))) - (4/3).*Ghm;
        Z = (Ghm./6) .* (9.*Khm + 8.*Ghm) ./ (Khm + 2.*Ghm);
        Gdry = 1./(((Phi/PhiC) ./ (Ghm+Z)) + (1 - Phi/PhiC)./(Gs+Z)) - Z;
    else
        % Phi >= PhiC
        Kdry = 1./(((1-Phi)./(1-PhiC)) ./ (Khm + (4/3).*Ghm) + (((Phi-PhiC) ./ (1-PhiC))./(Ks + (4/3).*Ghm))) - (4/3).*Ghm;
        Z = (Ghm./6) .* (9.*Khm + 8.*Ghm) ./ (Khm + 2.*Ghm);
        Gdry = 1./((((1-Phi)./(1-PhiC)) ./ (Ghm+Z)) + ((Phi-PhiC) ./ (1-PhiC))./(Gs+Z)) - Z;
    end
end

% Moduli saturated rock
num = (Phi.*Kdry) - ((1 - Phi).*(Kf .* Kdry))./Ks + Kf;
den = (1 - Phi).*Kf + Phi.*Ks - (Kf.*Kdry)./Ks;
Ksat = Ks .* (num./den); 
Gsat = Gdry;

% Velocities
Vp = ((Ksat + (4.*Gsat./3))./ RhoB).^0.5;
Vs = (Gsat ./ RhoB).^0.5;
% Bulk Density
Rho = RhoB;

%% PLOT
% D = [1:size(Vp,1)];
% figure;
% subplot(151);plot(Vp,D);grid on;title('Vp');axis ij;
% subplot(152);plot(Vs,D);grid on;title('Vs');axis ij;
% subplot(153);plot(RhoB,D);grid on;title('Rho');axis ij;
% subplot(154);plot(Phi,D);grid on;title('Phi');axis ij;
% subplot(155);plot(Sh,D);grid on;title('Sh');axis ij;
% 
% figure;
% subplot(131);scatter(RhoB,Vp,20,Sh,'filled');grid on;axis([1.2 2.6 1 5]);
% subplot(132);scatter(Vp,Vs,20,Sh,'filled');grid on;axis([1 5 0 3]);
% subplot(133);scatter(RhoB,Vs,20,Sh,'filled');grid on;axis([1.2 2.6 0 3]);
% 
% figure;
% scatter(Sh,Vp);grid on;


